Esta é uma pequena análise de regressão linear usando R. Meu objetivo aqui foi revisitar alguns conceitos, criar algumas funções e visualizar alguns dados referentes ao tema.

Vou usar o dataset Fish quem contém diferentes medidas corporais de peixes registrados em um mercado. Este foi retirado do kaggle e armazenado no meu dropbox, para posteridade. Todo código utilizado para gerar essa análise está disponível no repo do projeto, em meu github pessoal.

Antes, entretanto, de prosseguir, existe uma pergunta fundamental a ser respondida: porque e quando usar uma regressão linear? A grosso modo, ela pode ser utilizada para estimar um determinado valor (um valor esperado), com base em outro. Nesse contexto, assume-se que a relação entre estes dados é uma função linear.

1 Explorando o dataset e suas variáveis

O primeiro passo, como sempre, é se familiarizar com o dataset e carregar os pacotes que vou utilizar:

# pacotes utilizados pontualmente deveriam entrar como ::

library(rmarkdown)
library(tidyverse) 
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(ggExtra)
library(paletteer) # pacote com paletas de cores

paged_table( 
fish <- read_csv("https://www.dropbox.com/s/n45vml0ayoq0omx/fish.csv?dl=1")
)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   Species = col_character(),
##   Weight = col_double(),
##   Length1 = col_double(),
##   Length2 = col_double(),
##   Length3 = col_double(),
##   Height = col_double(),
##   Width = col_double()
## )

Os nomes das colunas não são lá muito intuitivos, por isso os descrevo abaixo:

  1. Species: nome da espécie mensurada.
  2. Weight: Peso do peixe em gramas.
  3. Length1: comprimento vertical em cm.
  4. Length2: comprimento diagonal em cm.
  5. Length3: comprimento cruzado em cm.
  6. Height: altura em cm.
  7. Width: largura diagonal em cm.
# convertendo a coluna species para fator
# nao compulsório, mas boa prática

fish <- fish %>% mutate(Species=as.factor(Species))

Meu objetivo aqui vai ser usar o peso como variável de resposta ou seja, quero usar alguma(s) das demais variáveis para prever qual vai ser o peso do peixe, o que, por sua vez, vai indicar (a grosso modo) seu porte.

Como pontapé inicial, vou plotar o peso em função de cada uma das demais variáveis, isso vai me dar uma visão geral dessas correlações:

Para plotar tudo de uma vez alterei um pouco a tabela e criei um grid:

fish %>% # alterando a disposição de colunas para a plotagem
  pivot_longer(cols=names(fish)[3:length(names(fish))],names_to = "explanatory") %>% 
  # plotagem
  ggplot(aes(y=Weight,x=value)) + geom_point() +   geom_smooth(method="lm",linetype=2) +
  theme_bw() +
  ggtitle("Peso vs. altura, medidas de comprimento e largura.") +
  ylab("Peso (g) - variável de resposta") +
  xlab("Variáveis explanatórias, em escala livre.") +
  facet_wrap(~explanatory, scales="free")
## `geom_smooth()` using formula 'y ~ x'

É possível perceber que a correlação do peso com as demais medidas é positiva em todos casos, pela inclinação das retas de regressão (tracejadas azuis), em alguns mais e outros menos. Entretanto, fica difícil entender a robustez dos modelos visualmente. Sendo assim, vou computar esses valores, na forma de R2, em uma tabela e ordenar de maneira decrescente:

# função para extrair o r2 dos modelos
rsquared_mod <- function(x) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = fish)
    )[["r.squared"]]
  )
}
# lista de variaveis explanatorias que vou cruzar com o peso
explanatory <- names(fish)[3:length(names(fish))]
# tabela com valores de r2 dos modelos em ordem decrescente
tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod))) %>% 
    arrange(desc(r_squared))

Essa tabela me dá os valores de \(R^2\) entre o peso e cada uma das outras variáveis independentemente (não sendo covariáveis) para cada modelo simples .

Seguindo por essa lógica, as três medidas de comprimento são as que melhor descrevem o comportamento do peso. Ainda, a natureza semelhante dessas medidas nos leva a crer que possivelmente existe uma colinearidade nessa relação.

Vou desenvolver isso um pouco mais a frente, mas por hora, como vou fazer uma regressão simples, vou usar a variável Length3 como explanatória, dado que ela resultou no maior valor de \(R^2\).

Abaixo, os coeficientes do modelo:

lm(Weight~Length3,data=fish)
## 
## Call:
## lm(formula = Weight ~ Length3, data = fish)
## 
## Coefficients:
## (Intercept)      Length3  
##     -490.40        28.46

Segundo esse modelo, existe um incremento de cerca de 28 gramas no peso de um peixe conforme a largura cruzada do peixe aumenta em uma unidade, ou seja, 1 centímetro.

Abaixo, vou plotar o gráfico da regressão, já atribuindo cores diferentes aos pontos referenes a cada espécies da amostra (estou levando todas em consideração):

2 A qualidade da regressão

fish %>% ggplot(aes(y=Weight,x=Length3)) + 
  geom_point(alpha=0.7,aes(color=Species)) + 
  geom_smooth(method="lm", linetype=2) +
  theme_bw() +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  labs(color='Espécie') +
  ggtitle(" Relação entre peso e comprimento cruzado,\n com cores referentes às espécies das observações.") +
  scale_color_paletteer_d("pals::alphabet") +
  theme(plot.title = element_text(size = 10))
## `geom_smooth()` using formula 'y ~ x'

Um fator interessante de se notar é a provável alta influência dos pontos associados à espécie Pike (em roxo), na inclinação e erro associado da reta, devido ao alto valor tanto de peso quanto de comprimento cruzado de alguns de seus pontos. Essa espécie provavelmente tem maiores valores, em média, de distância de cook para suas observações, o que vai influenciar no comportamento da reta. Podemos checar isso sumarizando o modelo e esses valores, por espécie:

mod <- lm(Weight~Length3,data=fish)

augment(mod,data=fish) %>% 
  group_by(Species) %>% summarise(mean_cooksd=mean(.cooksd)) %>% 
  arrange(desc(mean_cooksd))

A distância de Cook mede a influência das observações na inclinação e erro associado à reta de regressão, ou, seja, a influência destes nos valores ajustados: quanto maior esse valor, mais a inclinação da reta e seu erro está sendo definido por essas observações.

Como previsto, a espécie Pike de fato tem a maior média para os valores de distância de Cook associada a seus pontos, seguida pela espécie Smelt (verde escuro), que no gráfico tem uma concentração de pontos, relativamente alta, próximos de zero. Esses pontos, em ambas espécies, se distanciam bastante da reta de regressão, explicando o comportamento das médias de distância de Cook associadas a elas.

Vou checar mais a fundo a influência destas espécies na reta e no valor de \(R^2\).

Lembrando que o modelo tem \(R^2\) de 0.85, vou verificar qual valor que \(R^2\) assume retirando essas espécies dele:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.) %>% summary() %>% {.$r.squared}
## [1] 0.9109848

Retirando as espécies, a robustez do modelo aumenta significativamente. A inclinação da reta também se altera, conforme podemos ver na mudança do coeficiente associado à variável explanatória de comprimento cruzado:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.)
## 
## Call:
## lm(formula = Weight ~ Length3, data = .)
## 
## Coefficients:
## (Intercept)      Length3  
##     -656.48        34.14

Para enxergar essa alteração melhor ainda, podemos sobrepor as retas de cada modelo:

df_alt <- fish %>% filter(!(Species %in% c("Pike","Smelt")))
mod2 <- lm(formula = Weight ~ Length3, data = df_alt)


df_alt %>% ggplot(aes(y=Weight,x=Length3)) + geom_point(alpha=0) +
  geom_abline(intercept = -656.48, slope = 34.14, linetype=2,color="blue") +
  geom_abline(intercept = -490.40, slope = 28.46, linetype=2,color="red") +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  ggtitle(" Gráfico de pontos com reta de regressão (em azul) sem as espécies Pike e Smelt.\n A reta em vermelho é derivada do modelo com todas espécies.") +
  theme_bw() +
  theme(plot.title = element_text(size = 10))

Visualizando dessa forma fica mais claro o comportamento da nova reta em comparação com a reta associada ao primeiro modelo (tracejada vermelha). Essa inclinação indica uma correlação mais forte entre as variáveis na ausência das espécies Pike e Smelt, fato reforçado quando se calcula os índices de pearson para os dois casos:

tibble(df=c("df_original","df_sem_especies_ruido"),
       indice_pearson=c(
         cor(fish$Weight,fish$Length3),
         cor(df_alt$Weight,df_alt$Length3)
       ))

Finalmente, sem aprofundar muito, existe outra maneira prática de verificar o quanto os pontos das espécies podem estar influenciando na qualidade da regressão, que é plotando as curvas de regressão por espécie:

fish %>% ggplot(aes(y=Weight,x=Length3)) + geom_point() + geom_smooth(method="lm", linetype=2) +
  ylab("Peso (g)") +
  xlab("Comprimento cruzado (cm)") +
  coord_cartesian(ylim = c(0,1500)) + 
  ggtitle(" Peso vs. comprimento cruzado, por espécie da amostra.") +
  theme_bw() +
  theme(plot.title = element_text(size = 10)) +
  facet_wrap(~Species)
## `geom_smooth()` using formula 'y ~ x'

A espécie Pike é a única com valores de peso maiores que 1500. Já a espécie Smelt tem apenas valores próximos de zero, reforçando um pouco da influência destas na qualidade da regressão. Essa informação é redundante com o gráfico onde os pontos foram pintados de acordo com as espécies de referência, mas é uma maneira distinta interessante de enxergar esse comportamento.

Esse tipo de fator levanta a questão: quais espécies escolher para estabelecer um modelo geral? Uma maneira de fazer isso seria calculando o \(R^2\) associado ao modelo escolhido por espécie, mas também indicando o número de observações por espécie que existe nos dados:

# funcao para extrair os valores de R2 dos modelos
rsquared_mod <- function(x,data) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = data)
    )[["r.squared"]]
  )
}

# lista de vars explanatorias
explanatory <- names(fish)[3:length(names(fish))]

# tabela com valores de r2 dos modelos
tabela_r <- tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod,data=fish)))

# tabela ordenada

tabela_r %>% arrange(desc(r_squared))
temp <- function(sp){
  data <- fish %>% filter(Species==sp)
  return( 
  tabela_r <- tibble(Species=rep(sp,5),
                     variables=explanatory,
                     r_squared=unlist(map(explanatory,rsquared_mod,data=data)))
  )
}

temp2 <- map(unique(fish$Species),temp)

temp3 <- do.call("rbind",temp2) 

left_join( 
temp3 %>% filter(variables=="Length3") %>% arrange(desc(r_squared)),
fish %>% group_by(Species) %>% summarise(n=n()),
by="Species"
)

Quando analisadas em separado, a espécie Pike (que apresentava altos valores de distância cook, em média), não apresenta um \(R^2\) tão baixo. Mas é importante observar que o número de observações da maior parte das espécies da amostra é bem baixo (menor que 20). Apenas as espécies Perch e Bream tem mais de 20 observações. A espécie Whitefish, que contém o maior valor de \(R^2\) em seu modelo, contém apenas 6 observações. Se a ideia fossse utilizar um modelo generalista, seria razoável manter apenas espécies número alto de observações (20 talvez) ou aumentar o número de medidas das espécies subamostradas.

3 Colinearidade de variáveis

Como observado, todas medidas de comprimento apresentaram altos valores de \(R^2\) quando correlacionados com o peso. Recapitulando a tabela:

# funcao para extrair os valores de R2 dos modelos
rsquared_mod <- function(x) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = fish)
    )[["r.squared"]]
  )
}

# lista de vars explanatorias
explanatory <- names(fish)[3:length(names(fish))]

# tabela com valores de r2 dos modelos
tabela_r <- tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod)))

# tabela ordenada

tabela_r %>% arrange(desc(r_squared))

Esse comportamento pode ser um indicativo de colinearidade, o que faria sentido dada a natureza similar dessas medidas, em termos de proporção. Isso pode ser conferido computando os fatores de inflação de variância (Vif) para cada variável, em um modelo que leva em conta todas variáveis:

library(car)

fish %>% 
  select(-Species) %>%
  lm(Weight~.,data=.) %>%
  vif()
##    Length1    Length2    Length3     Height      Width 
## 1681.49649 2084.25783  422.46825   14.57009   12.27536

Os altos valores de VIF para as medidas de comprimento indicam que essas variáveis são, de fato, redundantes para o modelo, não fazendo diferença qual é utilizada: ambas tem boa perfomance para explicar os valores de peso. Explicando superficialmente, o Vif indica a severidade da colinearidade de variáveis.

Esse comportamento pode ser visualizado se plotarmos, por exemplo, o peso em função de dois comprimentos quaisquer:

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
temp <- fish

temp2 <- lm(Weight~Length3+Length1, data=temp)

p <- plot_ly(data = temp, z = ~Weight, x = ~Length3, y = ~Length1, 
             opacity = 0.6, colorbar = list(title = "Peso previsto")) %>%
  add_markers()

cf.mod <- coef(temp2)


x1.seq <- seq(min(temp$Length3),max(temp$Length3),length.out=1000)
x2.seq <- seq(min(temp$Length1),max(temp$Length1),length.out=1000)

z <- t(outer(x1.seq, x2.seq, function(x,y) cf.mod[1]+cf.mod[2]*x+cf.mod[3]*y))

p %>% # da para eu criar uma escala personalizada
  add_surface(x = ~x1.seq, y = ~x2.seq, z = ~z, 
              showscale = TRUE,colorscale="Viridis") %>% 
  layout(scene = list(xaxis = list(title = "Comprimento cruzado (cm)"), 
                      yaxis = list(title = "Comprimento vertical (cm)"),
                      zaxis = list(title = "Peso (g)")))
## Warning: 'scatter3d' objects don't have these attributes: 'colorbar'
## Valid attributes include:
## 'type', 'visible', 'showlegend', 'legendgroup', 'opacity', 'name', 'uid', 'ids', 'customdata', 'meta', 'hoverlabel', 'stream', 'transforms', 'uirevision', 'x', 'y', 'z', 'text', 'texttemplate', 'hovertext', 'hovertemplate', 'mode', 'surfaceaxis', 'surfacecolor', 'projection', 'connectgaps', 'line', 'marker', 'textposition', 'textfont', 'hoverinfo', 'error_x', 'error_y', 'error_z', 'xcalendar', 'ycalendar', 'zcalendar', 'scene', 'idssrc', 'customdatasrc', 'metasrc', 'xsrc', 'ysrc', 'zsrc', 'textsrc', 'texttemplatesrc', 'hovertextsrc', 'hovertemplatesrc', 'textpositionsrc', 'hoverinfosrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

Os pontos se dispersam muito pouco, o que indica sobreposição dos mesmos e sua redundância. Caso essa colinearidade não existisse, não seria esperada uma “quase reta” dentro de um espaço tridimensional.

4 Conclusões

No geral, a primeira conclusão que pode ser tirada da análise, é que o comprimento é a medida que melhor descreve o peso dos peixes, independente de qual é utilizado, uma vez que as medidas de comprimento são colineares. Além disso, é uma medida que por si parece ser um bom preditor do peso dos peixes.

Algumas observações, relativas a determinadas espécies de peixes, representam pontos aberrantes e comprometem a qualidade da regressão, como foi visto quando computando valores de distância de cook. Além disso, a quantidade de observações por espécie de peixe varia bastante, com diversas espécies com número baixo de observações, fator que deve ser levado em consideração na implementação do modelo.